rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_3b
#-figure_3c
#-figure_3a
#-figure_3g
#-figure_3d
#-figure_3e
#-figure_3f


library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library('gridExtra')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library('lfe')
library("ggrepel")
library("lemon")
library("sandwich")
library("multiwayvcov")
library("lmtest")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("geojsonsf")

#############################
#Function for calculating SE#
#############################

# Calculate robust confidence intervals
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))[, 2]
  return(res)}

pvalue_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))[, 4]
  return(res)}

ci_robust_cluster <- function(x){
  res_ci<-coefci(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))
  return(res_ci)}


data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")


#Step1. Turning some variables to numeric
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
data_cov_mob$spo_pop<-data_cov_mob$gspo/data_cov_mob$pop_total*1000
data_cov_mob$kul_pop<-data_cov_mob$gkul/data_cov_mob$pop_total*1000
data_cov_mob$nat_pop<-data_cov_mob$gnat/data_cov_mob$pop_total*1000


data_cov_mob$soz_pop<-data_cov_mob$gsoz/data_cov_mob$pop_total*1000
data_cov_mob$pol_pop<-data_cov_mob$gpol/data_cov_mob$pop_total*1000
data_cov_mob$inter_pop<-data_cov_mob$ginter/data_cov_mob$pop_total*1000
#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000


data<-data_cov_mob
data<-subset(data, week=="11")

m1<-lm(pct_afd~log(bridging_tot_pop+1), data=data)
summary(m1)

graph_brid_afd <- ggplot(data, aes(x=log(bridging_tot_pop+1), y=pct_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Bridging Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017", lim = c(0, 30))+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_brid_afd, file = "paper/graphs/figure_3b.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)

m1<-lm(pct_afd_2013~log(bridging_tot_pop+1), data=data)
summary(m1)

graph_brid_afd_2013 <- ggplot(data, aes(x=log(bridging_tot_pop+1), y=pct_afd_2013)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Bridging Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2013", lim = c(0, 30))+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_brid_afd_2013, file = "paper/graphs/figure_3c.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)


m1<-lm(pct_change_afd~log(bridging_tot_pop+1), data=data)
summary(m1)


graph_brid_change_afd <- ggplot(data, aes(x=log(bridging_tot_pop+1), y=pct_change_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Bridging Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017 - Pct. AfD, 2013", lim = c(0, 30))+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_brid_change_afd, file = "paper/graphs/figure_3a.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)


graph_bon_change_afd <- ggplot(data, aes(x=log(bonding_tot_pop+1), y=pct_change_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Bonding Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017 - Pct. AfD, 2013")+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_bon_change_afd, file = "paper/graphs/figure_3g.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)


graph_spo_afd_change <- ggplot(data, aes(x=log(spo_pop+1), y=pct_change_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Sport Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017 - Pct. AfD, 2013")+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_spo_afd_change, file = "paper/graphs/figure_3d.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)





graph_kul_change_afd <- ggplot(data, aes(x=log(kul_pop+1), y=pct_change_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Culture Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017 - Pct. AfD, 2013")+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_kul_change_afd, file = "paper/graphs/figure_3e.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)



graph_nat_change_afd <- ggplot(data, aes(x=log(nat_pop+1), y=pct_change_afd)) +
  geom_point()+
  geom_smooth(method=lm, se=FALSE, colour="black", size=2)+
  theme_bw() +
  scale_x_continuous(name = "Density of Nature Clubs / 1,000") +
  scale_y_continuous(name = "Pct. AfD, 2017 - Pct. AfD, 2013")+
  theme(axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16),
        axis.title=element_text(size=16),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(graph_nat_change_afd, file = "paper/graphs/figure_3f.jpg", 
       height = 15, width = 15, 
       units = "cm", dpi = 200)



